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Recent years have witnessed major advances in our understanding of nonequilibrium processes. 
The Jarzynski equality, for example, provides a link between equilibrium free energy differences 
and finite-time, nonequilibrium dynamics. We propose a generalization of this relation to non- 
Hamiltonian dynamics, relevant for active matter systems, continuous feedback, and computer sim¬ 
ulation. Surprisingly, this relation allows us to calculate the free energy difference between the 
desired initial and final equilibrium states using arbitrary dynamics. As a practical matter, this 
dissociation between the dynamics and the initial and final states promises to facilitate a range of 
techniques for free energy estimation in a single, universal expression. 


I. INTRODUCTION 

Free energy determination lies at the heart of nearly 
any application of statistical mechanics 0,i , the con¬ 
ventional methods being based on either the calculation 
of a partition function or the determination of work in 
a transition from one equilibrium state to another Q. 
In the latter case, the Helmholtz free energy Fo{A,P) 
of the initial equilibrium state with probability density 
oc exp —/ST-La is assumed to be known for a given 
Hamiltonian Ha with external parameters A and in¬ 
verse temperature ff = l/k's.T, where fee denotes the 
Boltzmann constant. Then, the free energy difference 
AFq = [Fq{B,I3) — Fo(H,/3)] corresponding to the tran¬ 
sition to another Hamiltonian Hs is estimated from the 
work as the external parameters are switched from A to 
B. To get an exact relation, the switching protocol is of¬ 
ten assumed to be either very fast, as in free energy per¬ 
turbation theory, or adiabatically slow, as in thermody¬ 
namic integration theory Q. A major breakthrough was 
achieved with the introduction of the Jarzynski equality 
(JE) 

g-/3AJ’o ^ (e-W”), (1) 

whereby the free energy difference AFq could be calcu¬ 
lated from the exponential average of the work W™ for 
any switching protocol of arbitrary speed. Here, the an¬ 
gular bracket denotes averaging over many repetitions of 
the switching protocol. The JE has lead to a plethora 
of new results in the context of nonequilibrium ther¬ 
modynamics and statistical mechanics |7Hl4l|. and there 
have been recent advances in the thermodynamics of con¬ 
trol [T5l - [l^ . predictio n self-replication [2l|, and in¬ 

formation processing |22h29I| . 

Before the introduction of the JE, Bochkov and Ku- 
zovlev had derived a similar relation (BKR) [jol - fs^ . 

1 = ( 2 ) 

where, surprisingly, is not equal to {W™ — AFq), as 
one might expect from the JE [Eq. ([T}]. (Note the def¬ 
initions of W™ and given below in Eqs. © and 


(IS|), respectively.) The apparent discrepancy between 
Eqs. [T] and [2] was resolved in [s^ by showing that 
these two relations correspond to two different conven¬ 
tions for defining internal energy, each leading to its own 
definition of work. In fact, by considering the dynamics 
under two different time-dependent conservative forces, 
/o and /i. Ref. [s^ presented a unified expression 

relating the free energy difference AFq to the different 
measures of work W™ and for the two forces /o and 
/i, respectively. The 0 subscript on AFq is meant to 
indicate that this free energy difference actually depends 
only on changes in /q; AFq is insensitive to changes in 
fi. We note that the detailed form of Eq. [3] differs from 
the corresponding formula presented in Ref. [s^, which 
assumes that the Hamiltonian is linear in fi. 

We point out here the surprising fact that Eq. © 
remains valid even if the system dynamics during the 
switching are not related to the two Hamiltonians Ha 
and Hb- Unlike the JE, where the dynamics during 
switching are derived from a time-dependent Hamilto¬ 
nian H{t) : Ha -T Hb, connecting the initial and fi¬ 
nal Hamiltonians, the combined JE and BKR formula 
[Eq. ([3])] implies that the intermediate Hamiltonian H{t) 
can be independent of them. More specifically, as long 
as the system is initiated in the equilibrium condition 
the dynamics during the switching can be governed 
by any modified Hamiltonian H'{t) = H{t) + Hi{t), for 
arbitrary Hi{t). To our knowledge, this aspect of the 
combined JE and BKR formula has not been appreci¬ 
ated before. 

Detailed studies have revealed the statistical quality 
of free energy estimation based on the JE [13, Isol - I^ . 
While slow driving protocols produce an effectively un¬ 
biased estimator, fast driving protocols induce far-from- 
equilibrium dynamics that often result in a bias. The 
convergence of a free energy estimator with respect to 
the number of independent samples is slow whenever the 
final phase space distribution of the system at the end 
of the protocol has a poor overlap with the final equilib¬ 
rium distribution Several strategies have been em- 
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ployed to improve the convergence of the JE estimator 
for fast driving protocols, including modifying the sys¬ 
tem dynamics to enhance the overlap between the ac¬ 
tual distribution of the system and and employing 
bidirectional protocols. However, while some previous 
studies have exploited specific forms of non-Hamiltonian 
dynamics in order to improve free energy estimation, it 
has not been clear what the optimal strategy should be, 
and a general framework unifying and extending previous 
results has been lacking. 

The main contribution of this paper is to introduce a 
generalization of the combined JE and BKR (Eq. (2), 
given by Eq. (EH) below, that is compatible with non- 
Hamiltonian dynamics during switching. In particu¬ 
lar, we extend the former strategy to arbitrary dy¬ 
namics, whereby the averaging in Eq. [T] is performed 
with the same initial equilibrium condition, p^, but 
the subsequent evolution is determined by potentially 
non-Hamiltonian dynamics, completely unrelated to the 
Hamiltonians T-Ca and Hb- 


II. ONE-DIMENSIONAL EXAMPLE 

Consider the isothermal dynamics of a particle of mass 
m constrained to move on a circle of circumference I. Let 
X and p denote its position and momentum, respectively, 
with the identification x -\-1 = x. There are two forces 
from the heat reservoir, damping {—jp/m) and noise {^t), 
the latter having the following statistical properties: 

(6) = 0, = (4) 

In addition, we consider two other forces: /o(a;;A) = 
—dxV{x; Ao), derived from some potential V{x; Aq) with 
external parameter Aq, and /i(x;Ai), with external pa¬ 
rameter Ai, which is not necessarily derivable from a po¬ 
tential. The dynamics of the particle are given by 

P 

X = —, 

m 

p = foix;Xo)fi{x;Xi) (5) 

m 


Even though our central result Eq. |37] is valid for ar¬ 
bitrary dynamics, for optimal estimation, the dynamics 
must be tailored such that the actual distribution of the 
system at the end of the dynamics is the same as the equi¬ 
librium distribution p^. Achieving this condition can 
make it possible to obtain accurate estimates after as few 
as one simulated transition. We derive the equation sat¬ 
isfied by such optimal modified dynamics and point out 
its relation to the so-called “escorted dynamics” [s^ ■ We 
emphasize that, whereas the optimal estimation strategy 
for a broad class of systems involves the use of escorted 
dynamics, it is typically not the case that one can com¬ 
pute the specific form these dynamics take for interesting 
non-equilibrium systems, pointing to the need for a more 
general framework such as ours. 

The manuscript is organized as follows. In Sec. H 
we illustrate our derivation of the modified JE compat¬ 
ible with arbitrary dynamics (Eq. 1371) for the simple 
case of the one-dimensional Langevin equation with only 
position-dependent forces. We first consider the under¬ 
damped case and then we describe some important sub¬ 
tleties of the overdamped limit associated with different 
stochastic integration schemes. In Sec. HI, we give a gen¬ 
eral proof of Eq. By construction, our general proof 
applies to situations involving many interacting Brown¬ 
ian particles. We show that, just like the JE, our result 
applies even when the dynamics take place in the absence 
of a thermal reservoir. In Sec. IV we emphasize how 
there is a clear separation between the dynamics and the 
end-point Hamiltonians Ha and in our new relation. 
In Sec. V, we derive an expression for optimal dynamics 
for free energy estimation based on Eq.[37]and discuss its 
relation to the so-called escorted dynamics. We conclude 
by comparing our result to some recent studies on work- 
fluctuation relations in the absence of detailed balance. 


where the dots over the variables x and p denote their 
time-derivatives. 

An equivalent way to describe the dynamics of the 
particle is via the Fokker-Planck equation for the phase- 
space probability density p(z,t), with z = {x,p): 

^=Cp=-V.-J, ( 6 ) 


where J = 


denotes the 


{^Jo + fi-lP-^dpj p 

phase space probability current. if the parameters 
{Ao,Ai} are held fixed in time, it can be shown that 
any initial distribution p(z, 0) relaxes to a unique sta¬ 
tionary distribution /9®(z) with Cp^ = 0 [s^. Note that, 
if the force /i is zero, the stationary distribution p® is the 
equilibrium distribution p®‘^(z; Xo,(3) oc exp with 

respect to the Hamiltonian Hxgiz) = p^ j{2m)-\-V{x\ Aq). 

Consider now a switching protocol specified by time- 
varying parameters {Ao(t), Ai(t)} with Ao(0) = A, 
Xo{t) = B, and arbitrary Ai(t). Following Refs. [nisi 
we now introduce two different notions of work, inclu¬ 
sive and exclusive. Inclusive work is applicable to only 
conservative forces while exclusive work is applicable to 
both conservative and nonconservative forces. Inclusive 
work done by the conservative force /o(a:; A) for a given 
protocol {Ao(t), Ai(t)} and over a trajectory {x{t),p{t)} 


IS 


W“ = / dt Ao(t) 


dV{x\Xo) 


9An 


( 7 ) 


x(t),Ao(i) 

The exclusive work done by the force fi{x] Ai) is 

= J dt ^ o/i(x(t); Ai(t)), 


( 8 ) 
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where the circle (o) on the right denotes Stratonovich 
multiplication [3^. According to the Feynman-Kac the¬ 
orem [IE [ 33 , the solution to the sink equation 

^ = Cg- hg, (9) 

with the initial condition g{z, 0) = Aq, /3) and arbi¬ 

trary phase space function /i(z, t), is given by the average 

g{z,t) = (^d{z{t)-z)e-fo'i*'Ht')^ (10) 

with h{t) = h{z{t),t) and S denoting the Dirac delta 
function. While Eq. OT is true of any h, if we consider 
the following particular form; 

h = pXodx,V + l3^h, (11) 

m 

its time-integral J dt h is equivalent to the sum W = 
{W^+Wf^), in units of 1//3, despite the fact that we have 
not used Stratonovich multiplication in the last term of 
Eq. (fTTl) (Appendix A). By direct substitution we can 
show that the unnormalized, time-dependent Boltzmann 
distribution 

is also a solution of Eq. m for the special choice of h given 
in Eq. [TTJ Combining Eqs. (fTUl) and (TT^ and integrating 
with respect to z at time t = t we get the following 
equation, 

g-/3AJ’o^^g-/3(w-+lvr)^^ (13) 

which is a special case of our more general result, 
Eq. (1^ . with AEh = Fo(i3,/?) — Fo(A,/3). Note that our 
approach — applying the Feynman-Kac theorem to the 
original protocol, as opposed to applying a Crooks-like 
fluctuation theorem to forward and reverse trajectories 
as an intermediate step — gives a much quicker 

derivation of Eq. dni) compared to previous approaches 
that have been brought to bear on this one-dimensional 
example HE- 


Q can also be thought of as the exclusive work done 
by the reservoir forces [IE- If we consider the follow¬ 
ing forms for the pi’s, 

pi(z) = p®‘i(z;B,/3), p 2 {z.) = p^^{z]A,j3), (16) 

then the fluctuation theorem ITU reduces to Eq. (fOl) . One 
just needs to use conservation of energy at the level of 
each trajectory: 

AHao = IF^ + VFr + Q, (17) 

derivable from the Langevin Eq. 


B. Subtlety in the overdamped limit 

Interestingly, our approach leads to a different integral 
fluctuation theorem than the entropy production fluctu¬ 
ation theorem in the overdamped limit, described by the 
following dynamics^: 

i = 7“^(/o + /i)+7“^6, (18) 

often useful in molecular simulations. In this limit, the 
two definitions of work W™ and bFf’' remain essentially 
unchanged; only {p/m) in Eq. ([5]) needs to be replaced by 
X. From the fluctuation theorem for entropy production, 
Eq. (I14L valid also in the overdamped limit, one can show 
the validity of Eq. (fOl) (Appendix B). However, using the 
Feynman-Kac approach described above, one can derive 
the following relation (Appendix C) 

^-PAFS ^ ^g-Jdt/ioD)^ ( 19 ) 

hoD = /3Ao ^AoF - ^ , (20) 

where, unlike the underdamped scenario, the quantity in 
the exponent on the right of Eq. [TO] is not equal to the 
sum PW = /3(W'™ -I- Wf^). I.e., Eqs. ITOl and [TO] are not 
the same. Note that the free energy Fg(Ao,/3) in this 
context is the configurational free energy 


A. Relation to fluctuation theorem for entropy 
production 


There is a close connection between Eq. (US and the 
integral fluctuation theorem for entropy production in 
the framework of stochastic thermodynamics 


piMillpwX 

P2(Z(0)) / 


= 1, 


(14) 


where pi and p 2 are any two normalized distributions 
and Q is the heat supplied to the system. 


Q = 


/ 



m 



m 


(15) 


^-l3FS{A,p) ^ z^{\o,P) = J ( 21 ) 

This duality of integral fluctuation theorems, observable 
only in the overdamped limit, has been reported be¬ 
fore [13 , but only in the context of the Bochkov-Kuzovlev 
relation, Eq. |21 In contrast, our treatment proves the 
existence of this duality in the broader context we con¬ 
sider here. Our preliminary investigation suggests that 
the duality stems from the effects of the transformation 


^ We have assumed uniform temperature and friction coefficient in 
the medium so that there is no anomalous contribution to the 
entropy production. 
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fi —>■ —fi on the probabilities of trajectories. We reserve 
the detailed investigation for a future study. 

We note that, from the perspective of numerical com¬ 
putation of AFq, it is advantageous to use Eq. (HH) in¬ 
stead of Eq. because the latter involves cancellation 
of large integrals, which is numerically costly. The re¬ 
quirement of the cancellation is easily seen via the Ito 
representation of the sum W = {W™ + W“): 


W = dt 


+ — ^/o/l + -^dxfl + fi + ft/l 


,(22) 

The fourth term on the right is nonnegative and its inte¬ 
gral grows with time. Convergence of the right hand side 
of Eq. (fT51) can be achieved only if this growing integral 
is canceled at each power. 


III. FLUCTUATION THEOREM VALID EOR 
GENERAL DYNAMICS 

We now derive a new fluctuation theorem that is valid 
even for general non-Hamiltonian dynamics (Eq. (1371) '). 
Consider a system described by the phase space coordi¬ 
nates z = {x, p} and the Hamiltonian 

^a„(z) = |^ + ^(x;Ao), (23) 

where p is the magnitude of momentum p and Aq = 
{Aqi, Ao 2 , . •.} is a set of external parameters. For a 
system with many particles, all of mass m, we have 
P = (PiiP 2 ,---) and X = (xi,X2,...) for particles 1, 
2, and so on. If we couple the system to a thermal reser¬ 
voir of inverse temperature /3, quite generally, we can 
write down the equation to motion to be of the following 
form^; 


P 

X = —, 

m 

P = -VxE(x;Ao)-r^+St, (24) 

m 

with (St) =0 and (SfS^) = 2{T/j3)5{t — t'). Here, T is a 
positive definite matrix denoting the damping coefficient 
matrix and S is the noise vector. The phase space dis¬ 
tribution p(z,t) evolves according to the Fokker-Planck 
equation 


^=4p = -Vz-J, (25) 

with J = [(p/to, -VxE(x; Ao) - Tp - (r//3)Vp)p]. As 
in the one-dimensional example, the asymptotic solution 


^ One can use different mass values rrn for different particles with- 
out changing the final result. 


is the Boltzmann distribution: 

p®'i(z;Ao,/3) = (26) 

e-I^FoiXo,P) ^ Zo(Ao,/3)= J dze-^'F-^°. (27) 

We now add an arbitrary phase space velocity vector 

vi = [fix(z; Ai),fip(z; Ai)] (28) 

to the dynamics of the system with external parameters 
Ai = {All, Ai 2 , ...}, leading to the following modified 
dynamics 



TO 


P — fo + fip ~ r-(29) 

TO 

where we have defined 

fo = -VxE(x;Ao). (30) 

Such additional phase space velocity vectors arise in 
many different contexts: (i) for velocity dependent feed¬ 
back control, with Vi = (0, —Tp/to) for some stable ma¬ 
trix r [^,1^; (ii) for self-propelled active particles, with 
vi = [0 ,F(p/to) 4- f(t)] for some odd function F and a 
generic function f [42l |: and (iii) for escorted, simulation 
dynamics, with Vi = AoiUi(z; Aq) for arbitrary, con¬ 
tinuous phase space vector fields ■ Note that 

vi can arise either from real physical forces, as in cases 
(i) and (ii) above, or from artificial dynamics intended to 
facilitate computer simulation and sampling of a system, 
as in case (iii). Note also that Vi does not generally fol¬ 
low from any Hamiltonian. The addition of Vi leads to 
a modified Fokker-Planck operator 

C = Co + Ci, £ip =-Vz • (vip), (31) 

leading to a modified stationary distribution p®(z), = 

0. Even when vi has a physical origin, i.e., it is of 
the form vi = [0, fip(z; Ai)] for some physical force 
fip(z; Ai), the stationary distribution p® may be un¬ 
known if the force is not derivable from a potential. 

Consider now initiating the system at the equilibrium 
distribution p®'^(z; Aq, /3) and driving the system accord¬ 
ing to some protocol A(t) = {Ao(t), Ai(t)}, as Aq varies 
from A to B. We wish to calculate the free energy differ¬ 
ence AFq = [Fo(B,/3) — Fo(A,/I)]. At any point along 
any trajectory z(t), the inclusive power by the original, 
conservative forces is given by 

= Ao-Vx„U(z(t);Ao), (32) 

and the exclusive power by the additional force 
fip(z(t); Ai) is given by 

fr/’'= fip(z(t); Ai) op(t)/TO. (33) 

However, unlike the one-dimensional case, the average of 
the exponential of minus the sum /3(IT™ -I- does 
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not give the free energy change A_Fo in this general case. 
We need to consider additional terms. To see this, let us 
begin with the unnormalized distribution 

By substituting it into the sink equation 

^ = Cg- hg, (35) 

and requiring that the right hand side of Eq. (l34l) is a 
solution of Eq. (1551) . we obtain the following expression 
for h: 

h = p + wr) - V, • Vi + • VxE. (36) 

Finally, applying the Feynman-Kac theorem to the sink 
equation (1551) with h as defined in Eq. (1551) . and combin¬ 
ing with Eq. (1541) . we get 

g-/3AFo ^ ^g-/3(W,)"-Hrvr)+/dt(V,-vi-/3fi,,.V,,V)^^ ^ 37 ^ 

which is our main result. 

Eq. (1571) generalizes Eq. (1T51) to the situation where the 
additional field Vi = (fix, fip) need not involve solely 
position-dependent forces. In particular, Eq. dsa in¬ 
cludes as a special case Brownian dynamics under elec¬ 
tromagnetic forces [ 4 ^. In this case, surprisingly, the 
integral in the exponent on the right of Eq. (157)) drops 
out leading to Eq. m, as observed in [d^. In a general 
scenario, both terms in the integral are non-vanishing. 
The first term, / dt Vz • Vi, accounts for the phase space 
contraction if the additional velocity term Vi is dissipa¬ 
tive. In the continuous feedback literature, this term has 
been referred to as entropic pumping [41| . The mean¬ 
ing of the second additional term, —f dtfix • VxE, is 
less transparent, probably because such terms do not ap¬ 
pear to have any physical origin, though they can arise 
in artificial, simulation dynamics [s^. 


Eq. 1361) over many repetitions of the protocol A(t) = 
{Ao(t), Ai(t)}. We can rewrite the integral (1//3) f dth 
along any phase space trajectory z(t) as 

= Jdt [(Ao • Vao + fip • Vp + fix • Vx) Hxo - (1//3)Vz ■ vi] 
= Jdt - (l//3)Vz • 

=7^b(z(t)) -'Ha(z(0)) - (1//3) J dtVz-z, (39) 

where we have used Eqs. [361 ESI and [55] in the first line 
and Eq. (55] in the second line. Note that, because the 
evolution of the system takes place in the absence of a 
thermal reservoir, the evolution is deterministic — if we 
know the initial phase space coordinate z(0), we know 
the future trajectory z(t > 0) for any given protocol 
{Ao(t), Ai(t)}. As a result, the integral (1//3) J dth can 
be treated as a function of just z(0). In particular, we 
can rewrite the average (exp(— f dth)) as 

J dz(0)p®‘i(z(0);A,/3)e-^‘^*^. (40) 

We can simplify Eg. l40l further: 


n p—/SW aIzIO ))-/dth 

(41a) 


(41b) 


(41c) 

^o(B;/3) 

^o(A,/3) 

(41d) 

g-/3[Fo(B,/3)-Fo(A./3)]^ 

(41e) 


A. Thermally isolated dynamics 

Just like the original JE, a feature of Eq. (1571) is that it 
is valid even when the dynamics during the switching are 
thermally isolated Q. In this case, the system is initiated 
in the same equilibrium distribution, p®‘i(z; A,/?), but 
the subsequent evolution does not involve the reservoir 
terms in Eq. [551 i.e, the system evolves according to the 
following dynamics: 


In line l41al we have rewritten Eq.HOl in line l41b[ we have 
used Eq. |55| for the integral f dt h; in line I41cl we have 
used the fact that the Jacobian |clz(r)/clz(0)| is given by 
exp (f dt Vz • z); and in the last two lines we have used 
the definitions in Eq. 1271 This completes our alternate 
derivation of Eq. |57|for thermally isolated dynamics. 


IV. ARBITRARY DYNAMICS 


X = — -f fix, p == fo + fip. (38) 

m 

The proof based on Feynman-Kac theorem still applies 
with a modified Fokker-Planck operator [s^- How¬ 
ever, the following derivation provides more insight. We 
evaluate the average (exp(—f dth)) (with h given by 


We want to emphasize that the dynamics during the 
protocol Aq (t) can be completely independent of Haq and 
we will still have Eq. (57] Consider an additional field Vi 
of the form 


fix = — 


m 


fip — 


5E(x;Ao) 

dx 


+ fip, 


(42) 
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for any given \o(t) and arbitrary Vi = (fix,fip)- The 
system then evolves according to 

X = fix, 

P = fip-r^ + St, (43) 

m 

without any term related to the Hamiltonian 
and yet we will still recover the free energy difference 
AFq from Eq. [33 (The reservoir terms E ^ and St will 
be missing in context of thermally isolated evolution of 
Sec. IIII Al l This indicates an interplay between the dy¬ 
namics and the quantity to average on the right hand side 
of Eq.|33which keeps the left hand side intact. This level 
of flexibility in choosing the dynamics seems not to have 
been appreciated before. Another benefit of the current 
approach is that we can quickly derive Eq. (1371) without 
going into detailed considerations of path integrals and 
conjugate processes [l^. Is^. l45j| . 


V. OPTIMAL DYNAMICS 

The dissociation between the dynamics and the initial 
and final equilibrium states promises to facilitate a range 
of techniques for free energy estimation in a single, uni¬ 
versal expression. Indeed such an instance has been seen 
before!^ for a special class of additional phase space ve¬ 
locity vector vi referred to as escorted dynamics. In the 
presence of a single time-dependent parameter Ao(t), the 
following form of vi was chosen, 


vi = Aou(z; Ao), (44) 

and it was shown that with an appropriate choice of 
u(z; Ao) it is possible to vastly improve the statistical 
quality of the free energy estimator based on Jarzynski- 
like relations. The essential idea behind choosing the 
appropriate dynamics was to ensure that the distribu¬ 
tion of the system under the modified dynamics, p(z;t), 
evolves close to the time-dependent equilibrium distribu¬ 
tion p®'i(z; Ao(t)). In fact, an exact equation was pro¬ 
posed for the optimal choice of u(z; Aq) by requiring that 
the time-dependent distribution is exactly the same as 
p®'^(z; Ao(t)). In this section, we consider the case of more 
than one time-dependent parameter. We see that the op¬ 
timal dynamics for free energy estimation can be recast 
in a generalized version of Eq. jUjin an extremely general 
setting. 

In order to obtain the optimal dynamics, such that 
a single instantiation is sufficient to yield an accurate 
estimate of the free energy difference, we can impose the 
condition that p®a( 2 ;Ao(t)) is a solution of the modified 
Eokker-Planck equation 


dt 


= -Cp, 


(45) 


(see Eos. [351 and 1311) then, after some algebra, we get the 


following equation^ for optimal Vi(z,t), denoted by v*: 

[Hxo - Ao(Ao;/3)]. 

In this case, we can obtain a solution of a form similar 
to Eq. |33| by considering an additional phase space ve¬ 
locity field v*j for each external parameter Aoi, be., by 
considering 

Ao), (47) 

i 

with each additional field u* now satisfying the following 
equation (under the assumption that all of the Aoi (f) are 
non-zero) 

V^-u*-Pn*-V,nx, = -Pdx,, [Uxo - Tl)(Ao;/3)] • (48) 

Equation |3S| can be simplified further due to the form 
of the Hamiltonian given in Eq. (33] Consider the nota¬ 
tion u* = (uj;,j,u*p). We can consistently assume that 
the momentum components, u*p, are zero and get the 
following equation for u*^(x; A): 

Vx ■ u*, - /3uL • VxH(x; Aq) = (E - Fq) . (49) 

As may be seen from Eqs. |46l |48l and |49l equations 
for the optimal dynamics are complicated and they even 
involve the free energy itself that we are trying to calcu¬ 
late. Clearly, it is extremely unlikely that one could de¬ 
rive the optimal dynamics in all but the simplest cases. 
Nonetheless, as for the case of escorted dynamics, the 
current approach provides insight into how to choose Vi 
such that free energy estimation is enhanced. The fact 
that escorted dynamics are already sufficiently powerful 
to provide the optimal dynamics for estimation might 
seem to suggest that there would be no practical benefit 
to developing tools compatible with a broader class of 
dynamics. However, note that we would already need to 
know the free energy difference, the very quantity we are 
trying to estimate, in order to solve Eq. (1351) in nearly any 
physical system of interest, which is why our more gen¬ 
eral relation Eq. (1371) promises to be useful for efficient 
free energy estimation. 


VI. DISCUSSION 


The current work should be contrasted with the stud¬ 
ies described in Refs. [4^.lisl lisl - fs^ . In only veloc¬ 
ity dependent additional forces were considered whereas 
in Refs. IH-lsS only position dependent forces were 
only those additional forces were 


_,|48| 
considered. In 


® As before, we are assuming isothermal dynamics. However, we 
might also consider time-dependent temperature by considering 
/9 as a parameter in the set Aq . 
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considered for which the steady state distribution was of 
the Boltzmann form, which is not the case in the cur¬ 
rent work. In [s^, [1^, the authors started with a generic 
Langevin equation, not derived from a Hamiltonian, and 
tried to build a thermodynamic theory for the dynamics. 
Their approach was based on the decomposition of their 
abstract dynamics into reversible and irreversible com¬ 
ponents. Our approach is complementary to theirs as 
we start from a given Hamiltonian and then add forces 
that may be nonconservative. Finally, in [ssj the au¬ 
thors propose a speed up of the calculation of free energy 
differences by utilizing the fact that violation of detailed 
balance can be used to accelerate the relaxation to steady 
states [ 5 ^. Given that Eq. ([37| is valid even in the ab¬ 
sence of detailed balance, it will be interesting to inves¬ 
tigate whether our new relation will lead to yet faster 
algorithms for calculating changes in free energy. 

Another interesting direction for future research con¬ 
cerns the further generalization of our results to include 
the determination of free energy profiles along reaction 
coordinates, as opposed to the free energy difference just 
between two given equilibrium states. Also, it remains 
to be seen to what extent the framework of bidirectional 
protocols developed for the Jarzynski relation and its 
generalization by Hummer and Szabo may be de¬ 
veloped for our new relation [s^, [13 ■ This could yield 
even more efficient approaches for calculating changes in 
free energy from limited data or simulations. 

In addition to the benefits of more efficient free en¬ 
ergy estimation techniques for physical systems that ac¬ 
tually obey Hamiltonian dynamics, the importance of 
developing a general framework compatible with non- 
Hamiltonian dynamics is highlighted by a recent example 
of an “active matter,” non-equilibrium system that 
cannot be handled by the JE, the combined JE and BKR, 
or any other previous generalizations that we are aware 
of. This particular system consists of a colloidal particle 
in contact with an “active bath” containing bacteria that 
themselves dissipate heat and create microscopic struc¬ 
ture in the solution, but other types of active matter 
systems have been observed [6J| . We hope and ex¬ 

pect that our framework will facilitate the study of active 
matter systems, as well as systems subject to continuous 
feedback, which to our knowledge are not correctly de¬ 
scribed by any previous generalization of the JE or the 
combined JE and BKR. 
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Appendix A: Irrelevance of Stratonovich scheme in 
underdamped dynamics 

Here we show the equivalence of the following stochas¬ 
tic integrals in the context of underdamped Langevin dy¬ 
namics: 

IF“ + wr = J dt (Ao^Ao ^ + /i o , (Al) 
^Jdth = J dt (AoaAoU + /i^) . (A2) 

The circle (o) on the right of Eq. m denotes 
Stratonovich multiplication |38l| - The Ito representation 
of integral (jATJ can be obtained through the following 
steps: 


J dt (Ao9aoV + /io£) 

fi{t) + fi{t + dt) p 


= Jdt 

= f dt 


^od\gV + 


(A3a) 


Ao^AqU + fi{t) -h 


p fiit + dt)-fi{t) p 


'.j dt (Ao5A„U + /i(t)£) 

[{dxfi)tdx + {^x^fl)^d\l] p 


dt ■ 


m 
(A3b) 


(A3c) 


In line IA3al we have used the definition of Stratonovich 
integration; in line lA3b1 we have simply rearranged terms; 
and in line IA3cl we have used the Taylor expansion of 
fi{t J- dt) to first order in dt; fi is a function of x and 
Al, both of which depend on time. We have neglected 
the higher order terms in the Taylor expansion for the 
reasons given below. 

The first integral in line IA3bl is equal to the integral 
in Eq. IA2I In the following, we argue that the other 
terms in line IA3cl are negligible. According to the under¬ 
damped Langevin equation [Eq. [S] of the main text], the 
differential dx varies as dt. If we assume the given pro¬ 
tocol Ai(t) to be continuous and smooth, the differential 
dAi also varies as dt. We may conclude that the whole 
integrand in the last integral of line [A3b| varies as dt, and 
therefore the integral is zero. The higher order terms in 
the Taylor expansion of line lA3a1 varv as higher powers of 
dt, and their contributions are zero as well. Gombining 
Eos. IA21 IA3al and lA3b[ we find that the two integrals in 
Eqs. lAll and IA2I are equivalent. 


Appendix B: Derivation of Eq. 1131 for 
one-dimensional overdamped dynamics 

Using the overdamped Langevin equation [Eq. [TH| of 
the main text], we can derive the following form of the 
first law of thermodynamics valid at the level of each 



























realization x(t): 


AE = Q + W^’^ + 


for an arbitrary function h(x, t) and initial condition 
g(x, 0) = exp[-^l/(a;, Ao)]/^'=(Ao(0),/3) is given by the 
(Bl) following expression 


where the energy E is equal to the potential energy V\ 
heat Q is given by Eq. [TS] of the main text, with p/m 
being replaced by i; and the two types of work W™ and 
are given by Eqs. [7]and[8]of the main text, respec¬ 
tively, again with the replacement p/m x. If we use 
Eg. IBlI in Eq. (TJof the main text, with pi and p 2 having 
the following forms. 


Pi oc e 


-PV{x-B) 


P2 oc e 


-PV{x-A) 

? 


(B2) 


g{x,t) = {5{x{t) - x)e (C3) 

with h{t) = h{x{t),t). Consider now the following ex¬ 
pression of h: 


h = hoD = ^ ( fofi + idj. 


(C4) 


after some cancellation of terms, we arrive at Eq. 1131 


Appendix C: Derivation of Eq. 1191 


For this particular choice of h, by direct substitution, we 
can show that the following expression of g also solves 
the sink equation IC2I 


The Fokker-Planck equation corresponding to the over¬ 
damped Langevin equation [Eq. |TS] in the main text] is 
given by 


— = Cod = —dxJoD, (Cl) 

where p{x, t) is the probability density at position x at 
time t, /loD is the overdamped Eokker-Planck operator, 
and Joo{x, t) = [7“^(/o + /i) - p is the prob¬ 

ability current density at position x at time t. (Other 
symbols have the same meaning as in the main text.) 
According to Feynman-Kac theorem, the solution of the 
sink equation 


dt 


CoDg - hg, 


(C2) 


gix,t) = 


,-PV{x-Xo) 


Z;/{A,P)= J 


Combining Eqs. IC3I and IC51 we get, 

p-PV{x-Xo) 


= ((5(x(t) -x)e-^o"di?^oD(t)). 


(C5) 


(C6) 


ZEiAP) 

Integrating both sides of this equation with respect to x 
at time t = r we get 


-PAFS ^ ^ , -/-dt/ioD(t)\ 

^5(A/3) ^ 


(C7) 


with AFq = F/j{B,P) — F§(B,P). Equation IC7I is the 
same as Eq. m of the main text. 
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